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For Int. J. Quantum Chem. 

Extending density functional theory (DFT) to an ab initio 
orbital functional theory (OFT) requires new methodology 
for nonlocal exchange and correlation potentials. This pa- 
per describes such modifications to a standard Dirac-Slater 
atomic program. Unrestricted Hartree-Fock (UHF) theory is 
extended by a modified Colle-Salvetti Ansatz for short-range 
electronic correlation. Results are reported for atoms He- 
Ne. Values of parameters needed for similar calculations on 
molecules and solids are reported. Implementation of nonlo- 
cal exchange and correlation for such extended systems, us- 
ing multiple scattering theory to connect independent calcu- 
lations in space-filling atomic cells, is discussed. 



I. INTRODUCTION 

Without invoking the electronic density as an interme- 
diate, variational equations for an independent-electron 
model can be derived as an orbital functional theory 
(OFT) [1], in agreement with Kohn-Sham equations in 
the local density approximation (LDA). In contrast to 
density functional theory (DFT) [2,3], an exact energy 
functional for OFT can be derived from many-body the- 
ory [4,5]. The orbital Euler-Lagrange equations (OEL) 
of OFT are determined by orbital functional derivatives, 
which take the general form of linear operators acting on 
occupied orbital functions of a model state. Because the 
exclusion principle requires independent normalization of 
the orbital partial densities, derivation of the OEL equa- 
tions requires orbital wave functions or densities to be 
varied independently [6]. For N electrons, there are N 
normalization constraints. Independent Lagrange multi- 
pliers (the orbital energy eigenvalues) are determined by 
these normalization constraints, as they are for Hartree- 
Fock and Kohn-Sham equations. 

In a density-based theory that incorporates the exclu- 
sion principle [7], Euler-Lagrange equations for the or- 
bital densities are determined by partial (Gateaux [8]) 
density functional derivatives, indexed by the orbital den- 
sities. Unless they are independent of the orbital index, 
these functional derivatives do not determine a unique to- 
tal (Frechet [8]) density functional derivative, equivalent 
to a multiplicative local potential function [6,9]. The re- 
quired Gateaux functional derivatives can be constructed 
from the corresponding orbital functional derivatives of 
OFT. The implied generalized Thomas-Fermi equations 



for orbital partial densities are operationally equivalent 
to the OEL equations [7]. 

Unless all orbital energies are equal [6,9], the OEL 
equations cannot be derived by a theory which considers 
only the total electronic density [10]. Equal orbital en- 
ergies conflict with the exclusion principle for the lowest 
state of any compact system with more than one elec- 
tron of each spin. This is easily verified in the example 
of the ls2s 3 S state of an atom with two noninteract- 
ing electrons of the same spin. Independent variation of 
orbital densities, required by the exclusion principle, de- 
fines Gateaux functional derivatives. It follows from this 
analysis that extension of Kohn-Sham theory beyond the 
LDA cannot be done in general without introducing non- 
local potentials [10]. The required nonlocal potentials are 
well-defined in OFT, formulated as an application of the 
many-body theory of electrons. 

The Schrodinger kinetic energy operator i = —5 V 2 
and the Fock exchange operator are well-known nonlo- 
cal "potentials" in standard theory. The present work 
includes nonlocal potentials for short-range correlation 
based on the correlation energy Ansatz of Colle and Sal- 
vetti [11]. The plan being followed is to implement these 
nonlocal potentials for calculations within atomic cells, 
then to use energy-linearized variational multiple scat- 
tering theory to match such local calculations together in 
large molecules and solids [12]. This methodology is out- 
lined here. Unrestricted Hartree-Fock (UHF) theory can 
be implemented for atoms or atomic cells by modifying 
any standard local-potential atomic Slater or Kohn-Sham 
program, to replace local exchange by the nonlocal Fock 
operator. A similar modification allows inclusion of the 
correlation energy functional considered here. Total and 
orbital energies computed with this nonlocal exchange- 
correlation model are reported for ground states of atoms 
He through Ne. 



II. NONLOCAL POTENTIALS IN ORBITAL 
FUNCTIONAL THEORY 

The simplest example of a nonlocal functional deriva- 
tive is provided by the Schrodinger kinetic energy orbital 
functional 



T = J2ni(i\i\i). 



(1) 



The notation here sums over occupied orbitals <j>i of a 
Slater-determinant reference state $, with occupation 
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numbers rij = 1 for i < N. Spin indices and sums are 
assumed, but suppressed in the notation. The total elec- 
tronic density of the model state p(r) — X)i n iP*( r ) ^ s a 
sum of orbital densities Pi(r) = </>|(r)</>j(r). The partial 
(Gateaux) density functional derivative 



ST . , <j>*A(j>i 



tiiSp. 



(2) 



is determined as an orbital-indexed local potential func- 
tion [7,10] by the orbital functional derivative 



ST 



t(f>i- 



(3) 



Applied to the functional V — n-i(i\v\i) , for an ex- 
ternal local potential v(r), this analysis determines the 
Gateaux derivative v v i(r) — v(r). Since this is inde- 
pendent of the orbital index, it reduces to a Frechet 
derivative, a multiplicative local potential. The implied 
generalized Thomas-Fermi equations for orbital densities 
i < N of noninteracting electrons are [7] 



v(t). 



(4) 



The Lagrange multipliers ej are to be determined such 
that J d 3 rpi = 1 for each orbital density. These equations 
are operationally equivalent to the noninteracting Kohn- 
Sham or OEL equations [7] 



(5) 



Two-electron functional are defined by U = E h + E x , 
where 

Eh = \ 'Y^ n i n 3{ i j\ u \ i j)'i E * = -\'Y^ninj{i3\v\ji), (6) 



and u — l/r\2- The Hartree functional Eh is an inte- 
gral of an explicit function of total density, which deter- 
mines the Coulomb potential Vh(r) as a Frechet deriva- 
tive. Similarly, if E xc is defined as in the LDA, the 
present analysis determines a Frechet functional deriva- 
tive, and verifies the LDA Kohn-Sham equations. The 
orbital functional derivative of the exchange functional 
E x defines the Fock exchange operator v x such that 
v x (f>i — —^j n j{j\ u \i) c t ) j- The functional derivative of 
U = Eh + E x defines u = Vh (r) + v x . Explicitly, 



SU 

niS(j)* 



^n 3 {j\u\j)<j) h 



(7) 



where u = u(l — P12), and V is the exchange operator 
for equal spin indices. 

Given (H — E)^f — for an N-electron eigenstate and 
any rule — > <& that determines a model state 3>, unsym- 
metric normalization ( < 1 > | V E') = ( < f ) | < E > ) = 1 implies E = 
{<f>\H\V) = E Q + E c . Here E = (<t>\H\<f>) =T + U + V 
is an explicit orbital functional, and E c = {<fr\HY& — $) 



defines the correlation energy. If Q — I — £? c — 

- $) = (<&\H\Q$). This implies an exact but 
implicit orbital functional [5] 

E c = -(<f>\H[Q(H -Eo-E c - i^Q^H^), (8) 

for 77 — > 0+. In practice, some parametrized approximate 
E c must be used, defining a correlation potential operator 
such that 



SE, 



(9) 



This defines an indexed local potential, the Gateaux 
derivative 



SE C 
n l Sp l 



v C i(r) 



(t>*v c <t>i 



(10) 



Defining universal functional F = E — V = T+ U + E c , 
and an operator T = t+u+i> c , the general OEL equations 
are 



SF 

m5(t>i 



= Tfa = - v(r)}<f>i, 



(11) 



reducing in the LDA to Kohn-Sham equations. The OEL 
equations imply generalized Thomas-Fermi equations [7] , 



= v f i(r) =€i- v(t). 



(12) 



If T is hermitian, the indexed potential v/i is the Gateaux 
derivative SF/mSpi. There is no implication in general 
that E x + E c defines a Frechet derivative. 



III. PLAN FOR LARGE MOLECULES AND 
SOLIDS 

The general case of nonlocal potentials determined ex- 
plicitly by the idempotent Dirac density matrix 



p(l,2)=^^(n)M!(ri) 



(13) 



is considered. In full-potential multiple scattering the- 
ory (MST) [13], local basis functions are constructed by 
integrating the Schrodinger or semirelativistic equation 
for specified orbital angular momentum I and energy 
e within the enclosing sphere (r = n) of each atomic 
cell. Following ideas of canonical energy-band theory 
[14,15], reviewed by Skriver [16], the energy-dependent 
radial wave function ug(e:r) is characterized by its loga- 
rithmic derivative Dg(e) = rsu' e (rs)/ue(rs), evaluated at 
the radius of a sphere whose volume equals that of a poly- 
hedral atomic cell. As discussed originally by Wigner and 
Seitz, this cellular wave function can continue smoothly 
across a cell interface if -C^(e) is negative, implying that 
Di(cb) = and Di(eA) = —00 define the lower and up- 
per energy limits of an energy band. The band center is 
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estimated to occur at Di{ec) = —I— I. An equivalent pa- 
rameter is pe (e) — D p^ ( ^, 1 , which varies nearly linearly 
over the width of a band. It is found that orbital wave 
functions are well-approximated by linear interpolation 
over energies es < e < ca- 

This behavior indicates that an initial self-consistent 
calculation should be carried out for spherically aver- 
aged potentials, subject to the orbital boundary condi- 
tion De(e) = —l—l at rs, in each inequivalent atomic cell. 
Because they satisfy a fixed boundary condition, the set 
of eigenfunctions u^{ec', r) can be extended to complete- 
ness within rs, defining a local basis set for the linearized 
variational cellular method (LVCM) [12]. The computed 
nonlocal potential must be extended out to the enclosing 
sphere r = ri, using the density matrix p constructed 
from the self-consistent cell orbitals, normalized to unity 
within rs- Orbital functions should be computed, using 
this fixed potential, for energies ec- Energy-derivative 
functions iii{cc',r) must be computed as a basis for ex- 
panding the LVCM global matching function [12,16]. 

At the UHF level of calculation, in each self-consistent 
iteration an average local exchange potential (Slater ex- 
change) is defined by 



(14) 



multiplying an antisymmetric model wave function. Us- 
ing coordinates q = — rj and r = j{ v i + r j) the CS 
Ansatz is 

£(r,q) = cxp(-/3V)[l - r(r)(l + l -q)]. (17) 
This ensures the limiting forms 



(1-01,-0 =r(r)(i + -?+•••) 

ZU^oo = 0. 



(18) 



In adapting this Ansatz to orbital functional theory 
(OFT) [1,5], it is desirable to retain the unsymmetric nor- 
malization condition ( < I ) | V I / ) = (3>|$) = 1, where is the 
correlated state, and $ is a reference state Slater determi- 
nant. Then the correlation energy is E c = ($>\H\$> — $) = 
J2i<j n i n jEij , a sum of electron-pair correlation energies. 
This suggests parametrization for each pair of occupied 



orbital functions 



such that 



E ij = -(ijHij(q)\ij). 
A parametrized form similar to CS is 



Z ij (q)=e X p(-fi j <?)[l->y ij (l + -q)]. 



(19) 



(20) 



an exact formula if v X i reduces to a Frechet derivative. 
This approximation is corrected by incremental inhomo- 
geneous terms 



{v x - v x }4>i 



(15) 



evaluated from the previous iteration. This procedure 
converts a local-potential algorithm into UHF, and is 
valid for the indexed correlation potential considered be- 
low. 



IV. SHORT-RANGE CORRELATION ENERGY 

Electronic correlation energy arises from two quite dif- 
ferent sources. At distances larger than an atomic radius, 
multipolar response produces correlation effects evident 
in polarization potentials and dispersion forces. Long- 
range correlation, not considered here, requires compu- 
tation of first-order multipole response pseudostates for 
each basis function ue(ec;r) [4]. A different approach is 
required for the short-range correlation due to the singu- 
larity of the intcrclcctronic Coulomb potential u = l/r\2- 
Expansion in r\i about such a singularity shows that an 
N-electron wave function must have specific cusp behav- 
ior [17,18] in order to cancel the singularity. The wave 
function must vary as 1 + \r\i + • • • . Colle and Salvetti 
(CS) [11] impose this cusp condition through a symmet- 
rical factor 



Parameter 7^ is determined by the normalization condi- 
tion ($|* - <f>) = 0, or {ij\^ij\ij) = for each pair: 



_ fa1cxp(-/3?V)|^') 
(y|exp(-^)(l+i g )|y)- 



(21) 



The free parameter can be chosen to minimize Bcthc- 
Goldstone (BG or IEPA) energy [19] for specified or 
can be treated as a semiempirical parameter and fitted 
to known correlation energies of atoms. 

This Ansatz can be incorporated into OFT by using 
Eq.(19) to define the model correlation energy. At any 
stage of the self-consistency iteration, parameter 7^ is de- 
termined by the consistency condition given above, and 
[3ij is either a fixed empirical parameter, or is to be up- 
dated by solving a 2-electron BG equation indexed by 
orthogonal occupied orbitals ij. The indexed local cor- 
relation potential (Gateaux functional derivative) is [10] 



where 



Vd(r) = 



#(r)&(r) ■ 



(22) 



(23) 



ni<j[l-£(ri,rj)], 



(16) 



a direct generalization of the operator u. Both Vh and 
v x are modified by short-range correlation, and antisym- 
metry is built in. Because terms j = i vanish, this 
Ansatz for Coulomb-cusp correlation does not produce 
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self-interaction. The indexed local potential v C i (r) is sin- 
gular at nodes of (j>i(r). As in Hartree-Fock methodol- 
ogy, v c 4>i should be treated as an inhomogeneous term in 
numerical solution of the orbital Euler-Lagrange (OEL) 
equations. 

Coupling of electron pairs is significant in standard CI 
methods for electronic correlation [19]. Although such 
coupling is inherent in the self-consistency of the pro- 
posed mean-field model, a more accurate extension of 
the method may be needed. A possible procedure is 
to implement the coupled electron-pair approximation 
(CEPA) [20] within each atomic cell, using the modified 
CS Ansatz as a closure formula for the CI expansion. 

A. Notes on integrals 

The normalized orbital basis functions are of the form 



with P-i(x) = 0,Pq(x) = 1. The elementary integrals 
required are of the form I n (/3; k) = J Q K e~P q q n dq. By a 
change of variables such that t — f3 2 q 2 , this reduces to 

= T %^ l)) A {n + l) ^ 2 ^ (30) 

in terms of the incomplete gamma function p(a, x) — 
1 — T(a,x)/T(a) [22], Sect. 6. 2. p(a,x) can be computed 
cffficiently using a power series for small x and a contin- 
ued fraction for large x. The continued fraction termi- 
nates if n is odd. 

V. CALCULATIONS ON LIGHT ATOMS 



<t> a (r)=N aX a(r)Y iama (6,cP), (24) 

where rx a { r ) is a numerical solution of the radial 
OEL equation. The normalization constants are N a — 
[J r 2 drXa( r )]~^ ■ Definite and indefinite integrals are re- 
quired for two-electron generalized potential functions 
F(q), where q 2 = r\ + r\ — 2rir 2 cos#. Integration over 
angles follows the standard derivation [21] of Condon and 
Shortley for F(q) = l/q. For normalized radial functions, 

(ab\F\cd) = c k (ac)c k (db)F k (ac; db), (25) 
k 

where c k denotes a Gaunt coefficient [21] and 

/>oo />oo 

F k (ac;db) = / r 2 dn / r\dr 2 
Jo Jo 

/ fe (^,r 2 )x:(ri)Xc(ri)x d (r 2 )x^(r2). (26) 
The factor f k is 

2k + 1 

f k (r u r 2 ) = ^— sm6d6P k {cosO)F(q). (27) 
1 Jo 

Given n and r 2 , cos 9 = and smQdQ = ^L. In 

internal coordinates, / fc (ri,r 2 ) is an integral of the form 

2fc + l r i+ - Pk{ rJ+A^ )qF{q)dq _ (28) 
2rir 2 J\ ri -r 2 \ 2rir 2 

These integrals are needed for qF(q) equal to l,q,q 2 
times the factor exp(— (3 2 q 2 ), so that the integrand of 
f k is this Gaussian factor times a polynomial in q. It can 
easily be verified that f k reduces to r</r> +1 for /3 — > 
if qF(q) = 1. The Legendre polynomial factors of the 
integrand are determined by the recurrence formula 

Ob J- 1 k 
Pk + i{x) = -^xPkW -j^Pu-xix), (29) 



For applications to molecules and solids, using multi- 
ple scattering theory and an atomic cell model, values of 
the parameters can be obtained by calculations on 
atoms. Results of such calculations, for light atoms He 
through Ne, are reported here. A numerical Dirac-Slater 
program [23] was modified as described above for UHF 
(exchange-only) and OFT calculations, the latter incor- 
porating the modified CS correlation energy functional 
described above. The program was used in its nonrel- 
ativistic mode. Angular coefficients in the total energy 
functional were computed such that effective potentials 
are spherically averaged, but retain a spin index. This 
equivalence restriction implies that radial orbital func- 
tions with indices £, m s are well-defined. 

For each atom considered, three sets of self-consistent 
calculations were carried out: for the He-like ion, for the 
Be-like ion (for N > 4), and for the neutral atom. To 
verify the computational method, computed UHF ener- 
gies are compared in Table(I) with established RHF en- 
ergies [24] and with total energies including correlation 
[25]. The OFT calculations were used to determine /3y 
parameters such that the computed total energies agreed 
with the "experimental" values shown in Table (I). Pair- 
indexed parameters /3i s i s were determined for He-like 
ions, parameters /3 2s2s for Be-like ions, and parameters 
/?2p2p for the neutral atoms. In each case, inner parame- 
ters were frozen and intershell parameters such as /3i s2s 
were scaled to the geometric mean of the corresponding 
intrashell parameters. Values of the latter that fit total 
energies of the ions and atoms considered are listed in 
Table (II). These parameters are the principal result of 
the present calculations, intended to define parametrized 
correlation functionals for extended systems. 

Self-consistent UHF orbital energies are tabulated in 
Table(III) and OFT orbital energies in Table(IV). Be- 
cause Janak's theorem [26] is valid in OFT, these ener- 
gies have a physical meaning in the context of a theory in 
which orbital occupation numbers are allowed to change 
continuously and to have fractional values. They are 
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derivatives of the total energy with respect to infinites- 
imal changes of the occupation numbers. Physical en- 
ergy differences, by implication, correspond to integrals 
of these derivatives over finite increments of occupation 
numbers. 



VI. CONCLUSIONS 

This paper considers an independent-electron model 
that incorporates a theoretically motivated Ansatz for 
correlation energy, expressed as a parametrized orbital 
functional. Computed results extend exchange-only the- 
ory (UHF) to a formalism parametrized by exact atomic 
ground-state energies. Parameters are obtained that 
make it possible to apply this formalism to calculations 
of the electronic structure of molecules and solids. 
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TABLE I. Total energies in Hartree units 



Atom 


RHF 


UHF 


6Xp 


He 


-2.8617 


-2.8617 




Li 


-7.4327 


-7.4328 


-7.4780 


Re 
DC 






1 A P.P.7A 


B 


-24.5291 


-24.5293 


-24.6541 


C 


-37.6886 


-37.6900 


-37.8466 


N 


-54.4009 


-54.4045 


-54.5890 





-74.8094 


-74.8136 


-75.0674 


F 


-99.4093 


-99.4108 


-99.7333 


Ne 


-128.5470 


-128.5470 


-128.9400 



TABLE II. Orbital (5 parameters 



Atom 


Is 


2s 


2p 


He 


0.83455 






Li 


1.42494 


1.16853 




Be 


2.01812 


0.52267 




B 


2.61268 


0.71557 


0.82167 


C 


3.20278 


0.89701 


1.08430 


N 


3.79877 


1.05841 


1.39511 


O 


4.39600 


1.20737 


1.55724 


F 


4.99306 


1.34754 


1.77434 


Ne 


5.57807 


1.47974 


2.00125 



TABLE III. UHF orbital energies (Hartree units) 
Atom ls a lsp 2s a 2sp 2p a 2pp 



He 

Li 

Be 

B 

C 

N 

O 

F 

Ne 



-0.91796 
-2.48668 
-4.73267 
-7.70036 
-11.34480 
-15.67067 
-20.70635 
-26.40567 
-32.77237 



-0.91796 
-2.46870 
-4.73267 
-7.68527 
-11.29972 
-15.58098 
-20.62807 
-26.35786 
-32.77237 



-0.19637 
-0.30927 
-0.54022 
-0.82450 
-1.16297 
-1.41447 
-1.66855 
-1.93040 



-0.30927 
-0.44175 
-0.57914 
-0.72580 
-1.07069 
-1.47217 
-1.93040 



-0.31671 
-0.43818 
-0.57092 
-0.67505 
-0.76666 
-0.85041 



-0.52144 
-0.67975 
-0.85041 



TABLE IV. OFT orbital energies (Hartree units) 
Atom ls a ls/3 2s a 2sp 2p a 2pfl 



He 

Li 

Be 

B 

C 

N 

O 

F 

Ne 



-0.94667 
-2.52396 
-4.81661 
-7.80392 
-11.46803 
-15.80892 
-20.86288 
-26.58113 
-32.96943 



-0.94667 
-2.50794 
-4.81661 
-7.79277 
-11.43171 
-15.73252 
-20.79561 
-26.53950 
-32.96943 



-0.19546 
-0.31196 
-0.53861 
-0.81782 
-1.15024 
-1.39967 
-1.65541 
-1.92025 



-0.31106 
-0.45076 
-0.59356 
-0.74332 
-1.08080 
-1.47252 
-1.92025 



-0.32004 
-0.44126 
-0.57450 
-0.67795 
-0.77354 
-0.86284 



-0.54422 
-0.69690 
-0.86284 



G 



